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Abstract 

We perform an analytical and numerical study of the crossover from the Joseph- 
son effect to the bulk superconducting flow for two identical one-dimensional su- 
perconductors, co-existing with a layer of normal material. 

A generalized Ginzburg-Landau (GL) model, proposed by S.J. Chapman, Q. 
Du and M.D. Gunzburger [1] was used in modeling the whole structure. When 
the thickness of the normal layer is very small, the introduction of three effective 
5-function potentials of specified strength leads to an exact analytical solution of 
the modified stationary GL equation. 

The resulting current density-phase offset relation is analyzed numerically. We 
show that the critical Josephson current density j c corresponds to a bifurcation of 
the solutions of the nonlinear boundary value problem coupled with the modified 
GL-equation. The influence of the second term in the Fourier-decomposition of 
the supercurrent density-phase relation is also investigated. 

We derive also a simple analytical formula for the critical Josephson current. 

1 Introduction 

If two superconductors are weakly coupled and have a phase difference A0 that is not or 
7r, a zero- voltage supercurrent flows from one to the other usually at a rate proportional 
to sin A(f), as predicted by Josephson [2] in 1962 for the case of junction consisting of 
an insulating oxide layer between two identical superconductors. A DC supercurrent 
can flow through such junction in the absence of a voltage difference between both 
superconductors in such a way that 

j a (A<f>) =j c sinA0. (1) 

After the discovery of the Josephson effect, it became clear that, apart from an in- 
sulating tunnel structure (SIS), any sufficiently short localized weak link such as a very 
short constriction in the cross-section of a superconductor, a point contact between two 
superconductors, as well as two superconductors separated by a thin layer of normal 
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metal, could be used as a Josephson junction, obeying the current-phase relations, usu- 
ally different from (f). This fact forced Licharev [3] and Waldram [4] to dwell upon 
a broad minded definition: a weak link is supposed to show a Josephson behavior if 
the supercurrent-phase relation is a single-valued and analytical function, which can be 
represented as a Fourier series 

oo 

i s ( A0) = J2 a n sin (nA0) . (2) 

n=l 

The crossover between an ideal Josephson behavior and an uniform superconducting 
flow was studied by solving exactly the usual Ginzburg-Landau (GL) equation for a 1-D 
superconductor in the presence of an effective 5-function potential of arbitrary strength 
(see, for example [5]). Recently, a modified GL type model has been formulated [1]. This 
model could equally well be applied to a boundary between different superconductors, 
superconductor-insulator, superconductor-normal metal. The purpose of our paper is to 
apply this modified GL model for calculating the supercurrent-phase relation and the 
crossover between Josephson behavior and uniform superconducting flow. 



2 Statement of the Problem 



The bifurcation analysis of superconducting solutions from the normal solutions for the 
1-D GL equations in the presence of an external magnetic field was considered in [6]. 
Here, we also consider the 1-D case, but on the basis of the generalized (modified) GL 
model [1]. 

By neglecting the external and the self-induced magnetic fields we have in the super- 
conducting domain x G (— oo, —d/2) and x G (d/2, oo) 
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and the boundary conditions at the N-S interfaces x = ± d/2 are 
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so that ] s = j n = j. Here [/] denotes, as usually [1], the jump of the enclosed function 
f(x) across the points x = ± d/2. The quantities e s and e n are equal to the charge of 
the superconducting charge carriers, and in the following are both equal to twice the 
electron charge. We are free to arbitrary choose one of the masses m s and m n . Usually 



m s is chosen to be twice the electron mass. This leaves parameter depending 

on the normal material. 

We suppose that a s = —\a s \, a n > 0, b s > 0, b n > and we define the coherence 
length 

y/2m s \a s \ ' 

as well as the dimensionless distances z — x/£, d — d/£, the order parameter ip(z) = 
t/j(z^)^b s /\a s \ : and the current density 



^|O s |O s |6 s 



In order to make further comparisons with other papers we introduce an equivalent 
representation of equation (6) 



where A = m s /nQe 2 s = /^oA|, and n = y/\a s \/b s is the equilibrium concentration of 
superelectrons, is the well-known London penetration depth, O — h/2e = h/e s is the 
magnetic flux quantum. 

With the definitions given in (5)- (6) our problem is stated as follows 

(1- = 0, (-00 < z < -d/2) U (d/2 < z < oo), (7) 

^"--^LV-^HV = 0, \z\<d/2. (8) 
m n \a s \ b s 

J=-):—[PV-W>~\, 
2m n 

3 Analytical solution for a thin normal layer 

Let us introduce the parameters 

Tfl n Tfl n (X n (X n Tfl n b n b n 

m = — , a = — j — 7 = m — , p = — = m— . 
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In both Eqns. (7), (8) we set ip(z) = R(z) exp [i<^(z)] and we find 
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R" + R-R 3 - — = 0, \z\ > d/2, (10) 



7 2 

R" - aR - (3R 3 - m 2 — = 0, \z\ < d/2. (11) 

It is not surprising that the case a — — 1, (5 — m — 1 corresponds to an uniform 
superconductor occupying the whole space (— oo, oo). 
Let us introduce the function 

1, \z\ > d/2; 



5{z; 1 -c) = 1 + (1 -c) 
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Here, 8(z) is the Heaviside function. Then Eqns. (10), (11) can be written in the 
whole space (— oo < z < oo) as 

J 2 

R" + 5 (z; 1 + a) R - 5 (z; 1 - (3) R 3 - 5 (z; 1 - m 2 ) — = 0. (12) 

If the thickness d — > +0, we have 9 {z — | ) — 9 {z + | ) = — dS (z) , so that the case 
of small d can be formulated as follows 

R" + [1 - 9l 8(z)} R-[l- g 2 5 {z)\ R 3 - [1 - g 3 5 {z)\ ^ = 0, (13) 

where we merely substitute g± = (1 + a) d, g 2 = (1 — (3) d, g 3 = (1 — m 2 ) d 
The solution of Eqn. (13 ) is found to be 

R 2 (z) = a + bt&nh 2 [u(\z\+ z )] (14) 



where a(2 - a) 2 = 8J 2 , < a < 2/3, 6=1- 3a/2, w = y^6/2, b = aB 2 , and the quantity 
2/o = tanh(-uzo) satisfies the following equation (0 < y < 1) 



v / 26S 2 y (l - y 2 ) = 9i(l + 5 2 y 2 ) " Ml + 5 2 y 2 ) 2 - a2(1 ^ 2l/g) - ( 15 ) 
If (yr 2 = and 33 = we recover Eqn. (15) from [5]. 

Now we will introduce the phase offset Aip. Eqn. (9) can be rewritten as 
so that for small d we have 

v / = ^[l-d(l-m)S(z)] (16) 

Due to the fact that the boundary conditions for the order parameter at z — > ±oo 
must be J A 

J R'(±oo) = 0, ip(z) = —z± —, z ^ ±oo, 

we derive from Eqn. (16) 

Jd(l-m) T f°°, 1 1 . , 

Here i?oo = R(±oo) > and R 2 (0) = a + by 2 = a(l + B 2 yl). Then by calculating the 
integral in Eqn. (17), we finally find 

Aip = 2{arctan5 - arctan( J By )} - jflv ( 18 ) 

a(l + B y ) 

If m = 1 this result formally coincides with Eqn. (16) in [5]. 



4 Numerical Results 



The Generalized Continuous Analogue of Newton's Method (see the survey by Puzynin 
et al. [7]) for solving the nonlinear ODE (12) on the finite interval z G (—L,L) with 
zero Neumann's conditions at the boundaries z = ±L and appropriate conditions at 
z = ±d/2, is applied. At each iteration the corresponding linear boundary value problem 
is solved numerically using the Finite Elements Method on a nonuniform grid, condensed 
to the boundaries z = ±d/2 of the layer. 

We note that Eqns. (9)-(ll), the boundary conditions at z = ±L, as well as the Veier- 
strass conditions at the points z = ±d/2, can be interpreted as a necessary extremum 
conditions for the free energy functional 

L 

F[R,(p] = J A(x, R, R', (p')dx + J [p(-L) - ip(L)} . (19) 

-L 

Here, the energy density A is given by 

if r' 2 + it>y 2 - r 2 + \r\ z i (-f , i) 

2 I £ ( R>2 + ^V 2 + aR 2 + I R A ) ,ze(-H) 

All numerical results from now on were obtained for L = 16 and width of the layer d = 
0.2. Therefore, the two superconducting layers conform to the finite intervals (-L, —d/2) 
and (d/2,L). The graphics displayed in Fig. 1 correspond to J(A<p) curves for four 



Figure 1: Some of curves J(A(p) for g 2 = and g 3 = 0. 

different values of gi (gi — 0, g\ — l,gi = 5, and gi = 10) at g 2 — and g 3 = 0. If 
the quantity g\ — (the corresponding curve is marked by □) the maximum is achieved 
at jdep = 2/3^/3 ~ 0.385 (uniform superconductor). For large g\ (gi — 5;gi — 10) we 
found results close to the ideal Josephson relation J = j c sin Aip, which will be analyzed 
more strictly below. We note that the numerical results, displayed on Fig. 1 are in good 
agreement with Fig. 2 in [5]. For each curve on this figure we denote jc = max J (Aip) 
when Aip/ir e [—1, 1]. 

For a given value of the current density J we found numerically two solutions, whose 
amplitudes R(z) and phases (p(z) are demonstrated on Fig. 2. For the first ("upper") 



Figure 2: The solutions for J = 0.1, g\ — 1, g 2 = and g 3 = 0. 



solution (marked by V) we have the phase offset Atp/ir G [—1, —j c ) and Atp/ir G (j c , 1], 
while for the second solution (marked by A) we have (Atp/ir G (—j c ,jc))- The first 
solution originates from the "uniform" solution R(z) = 1, tp(z) = 0, existing in the case 
when gi — 0, g 2 — 0, (? 3 = 0, and J = 0. 

The dependence of the free energy F(J) on the current density J for these two 
solutions is represented graphically on Fig. 3. This is a typical bifurcation diagram: in 
the point B we have J = j c , the two curves coalesce and acquire a common cusp. 



Figure 3: The critical current j c corresponds to a bifurcation point. 

Figures 4 and 5 are numerical investigations of the influence of the parameters g 2 ^ 
and g 3 ^ 0, respectively, on the J(Atp) curve. A comparison between Fig. 4 and Fig. 
1 at gi = 1 clearly shows that the influence of the parameter g 2 G [0, 1) on the current 
density is not very pronounced (a few percent), whereas a comparison between Fig. 5 
and Fig. 1 at g± — 1 proves that the variation of the parameter g 3 between and -3.5 
leads to a significant reduction of the maximum current density (approximately twice). 

These quantitative conclusions can be coupled with the Fourier decomposition of 
J(Atp) curves as given by Eqn. (2). We restrict ourselves only to the analysis of the ratio 
a 2 /ai of the first two Fourier coefficients. When 02/01 C 1 we have an approximately 
pronounced Josephson behaviour J ~ j c sin Atp = a x sin Atp. 



Figure 4: The Influence of the Parameter g 2 . 



Figure 5: The Influence of the Parameter g 3 . 



The ratio 02/01 as a function of the parameter gi is shown on Fig. 6. It is seen that 
for large values of the parameter g 1 (g 1 > 8) when the parameter g 3 = 0, the coefficient 
a 2 is less than 5% of a±. On the contrary, for small values of g\ we have a substantial 
weight of higher harmonics (for example, if g\ — 1 then the ratio a 2 /oi ~ 0.23). 

As it can be expected (see the curve marked by A), the influence of the alteration of 
the parameter g 2 on the Fourier coefficients is essential for small enough values only of 
the parameter g\ (gi < 4). On the other hand, taking into account the coefficient g 3 7^ 
(m^l) leads to a significant increase of the second term in Eqn. (1) (the corresponding 
curve 02/01 is marked by V). 

Let us consider the special case of small current J — > +0. In this case, from the usual 
sinusoidal Josephson relation (1) in linear approximation we get 



Eqn. (15) is simplified considerably if J = and reduces to the following equation 



where Y = y (J = 0). For small J we have b ~ 1, q = 2 J 2 , B = l/(y/2J), 1 > B' 1 and 



A(p = arcsin(J/j c ) ~ J/j c . 



(20) 



V2(l - Y 2 ) = 9l Y - g 2 Y 3 , 



(21) 



Figure 6: The Ratio a 2 /ai of the first two Fourier coefficients as a function of g±. 



the right-hand-side of Eqn. (18) can be replaced by a term, proportional to J 



A(p = 2J 
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(22) 



By coupling Eqns. (20) and (22) we derive an approximate estimation for the critical 
current 



(23) 



where Y is the smallest root of Eqn. (21). In the special case g 2 = 0, m = l(g 3 = 0), 
gi 3> \/2, Y Q = y/2/gi <C 1, we recover the result from [5] 

■ I^L _!_ 

Figure 7 shows the comparison between numerically obtained and theoretically cal- 



Fi gure 7: The Critical Current j c as a function of the parameter g x . 

culated curves j c (gi) by means of formula (23) for different values of parameters. We 
emphasize the agreement between the theoretical and numerically obtained relations. 



Concluding Remarks 



The physics of Josephson junctions is based on a usual sinusoidal supercurrent-phase 
difference relation. In the present paper, we show that, by taking into account different 
nonlinear terms in the normal and superconducting regions, many harmonics exist and 
the dependence J(A(p) of the current as a function of the phase offset is not sinusoidal. 
This effect follows from the numerical investigation and is seen also from the approximate 
analytical expressions given above. 

We show that the maximum current density j c represents a bifurcation point for the 
amplitudes of the supercurrent flow by a change of the current density J. 
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